arXiv: 1509.04629v2 [physics.comp-ph] 21 Nov 2015 


A minimization principle for the description of 
time-dependent modes associated with transient instabilities 

Hessam Babaee^’^, Themistoklis P. Sapsis^ * 

^Department of Mechanical Engineering, MIT 
^Sea Grant College Program, MIT 


Abstract 

We introduce a minimization formulation for the determination of a finite-dimensional, 
time-dependent, orthonormal basis that captures directions of the phase space associated 
with transient instabilities. While these instabilities have finite lifetime they can play a crucial 
role by either altering the system dynamics through the activation of other instabilities, or 
by creating sudden nonlinear energy transfers that lead to extreme responses. However, 
their essentially transient character makes their description a particularly challenging task. 
We develop a minimization framework that focuses on the optimal approximation of the 
system dynamics in the neighborhood of the system state. This minimization formulation 
results in differential equations that evolve a time-dependent basis so that it optimally 
approximates the most unstable directions. We demonstrate the capability of the method 
for two families of problems: i) linear systems including the advection-diffusion operator in a 
strongly non-normal regime as well as the Orr-Sommerfeld/Squire operator, and ii) nonlinear 
problems including a low-dimensional system with transient instabilities and the vertical 
jet in crossflow. We demonstrate that the time-dependent subspace captures the strongly 
transient non-normal energy growth (in the short time regime), while for longer times the 
modes capture the expected asymptotic behavior. 


1 Introduction 

A broad range of complex systems in nature and technology are characterized by the presence of 
strongly transient dynamical features associated with finite-time instabilities. Examples include 
turbulent flows in engineering systems (e.g. Kolmogorov [15] and unstable plane Couette flow 
[26] . reactive flows in combustion isam]), turbulent flows in geophysical systems (e.g. climate 
dynamics isnEni, cloud process in tropical atmospheric convection IMIES]), nonlinear waves 
(e.g. in optics ismi or water waves [33 EH uni mi), and mechanical systems [3a|5l[29l|35l[25]. 

These systems are characterized by very high dimensional attractors, intense nonlinear energy 
transfers between modes and broad spectra. Despite their complexity, the transient features of 
these dynamical systems are often associated with low-dimensional structures, i.e. a small number 
of modes, whose strongly time-dependent character, however, makes it particularly challenging 
to describe with the classical notion of time-independent modes. This is because these modes, 
despite their connection with intense energy transfers and transient dynamics, often have low 
energy and hence they are “buried” in the complex background of modes that are not associated 
with intense growth or decay but only with important energy. These transient modes often act 
as “triggers” or “precursors” of higher energy phenomena or instabilities and a thorough analysis 
of their properties can have important impact for i) the understanding of the system dynamics 
and in particular the mechanisms associated with transient features (see e.g. isniEiisi]), ii) the 
prediction and quantification of upcoming instabilities that are triggered through these low-energy 
dynamical processes (see e.g. [53l[T8l[T9|), hi) the control and suppression of these instabilities 
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by suitably focusing the control efforts in the low energy regime of these transient phenomena 
(seee.g. gOl ESI [IT] ). 

Transient dynamics is central in understanding a wide range of fluid mechanics problems. In 
the context of hydrodynamic stability, non-normality of the linearized Navier-Stokes operator 
can cause significant transient energy growth [inmii. It has been well-established now that 
the eigenvalue analysis fails to predict the short-time evolution of perturbations for convective 
flows [59]. Instead, the transient energy growth can be better understood by analyzing the 
pseudospectra of the linearized operator [58l[42]. Transition from laminar flow to turbulence, 
is another active area of research in fluid mechanics, to which, understanding the transient 
dynamical features is crucial. For a finite-amplitude disturbance, bypass transition has been 
observed in wall-bounded shear flows |49l [43] , in which case the non-normal growth of localized 
disturbances leads to small turbulent spots, bypassing the secondary instability process [27] . 
Recent computational and experimental studies also demonstrate the sudden transition from 
laminar to turbulent motion in pipe flows, where turbulence forms from localized patches called 
puffs [36l|6]. On the other hand, intermittent behavior is the hallmark of turbulent fluid flows. 
Turbulent ehaotie bursts appearing in spatially and temporally localized events can dramatically 
change the dynamics of the system [21]. Prototype systems that mimic these properties were 
introduced and analyzed in [32] [131 ESI [34] . Transient dynamics have a fundamental role in the 
intermittent behavior of passive tracers as well, where even elementary models without positive 
Lyapunov exponents [TT] [12] have been able to reproduce intermittent behavior observed in 
complex models. For such systems, the fundamental role of the random resonance between Fourier 
modes of the turbulent velocity field and the passive tracer has been recently illustrated m- 

In the context of uncertainty quantification and prediction a new family of stochastic methods 
relying on the so-called Dynamical Orthogonality condition was recently developed to deal 
with the strongly transient features of stochastic systems. The Dynamically Orthogonal (DO) 
Field Equations [46j [47] and Dynamically Bi-orthogonal (BO) equations [16] evolve a subspace 
according to the system stochastic PDF and the current statistical state of the system. Despite 
their success in resolving low-dimensional stochastic attractors for PDFS [44] [45] [48], these 
methodologies are often too expensive to implement for high dimensional systems (e.g. DO or 
BO require the simultaneous solution of many PDFs). In addition, in order to obtain an accurate 
description of the time-dependent dynamics many modes should be included in the analysis and 
the computational cost increases very rapidly (especially for systems with high complexity). 

Our aim in this work is to develop a method that will generate adaptively a time-dependent 
basis that will capture strongly transient phenomena. This approach will rely on system observ¬ 
ables obtained either through high fidelity numerical solvers or measurements, as well as the 
linearized equations of the system. The core of our approach is a minimization principle that will 
seek to minimize the distance between the local vector field of the system, constrained over the 
direction of the time-dependent modes, and the rate of change of the time-dependent modes. A 
direct minimization of the defined functional will result in evolution equations for the optimally 
time-dependent (OTD) basis elements. For systems characterized by transient responses these 
modes will adapt according to the (independently) computed or measured system history in a 
continuous way, capturing at each time the transiently most unstable directions of the system. 
For sufficiently long times where the system reaches an equilibrium we prove that the developed 
equations provide the most unstable directions of the system in the asymptotic limit. 

We demonstrate the developed approach over a series of applications, including linear and 
nonlinear systems. As a first example we consider the advection-diffusion operator where we show 
how the OTD basis captures the directions associated with the non-normal behavior. The second 
example involves the Orr-Sommerfeld/Squire operator that governs the evolution of infinitesimal 
disturbances in parallel viscous flows. Our goal here is also the computation of time-dependent 
modes that explain the transient growth of energy due to non-normal dynamics. The third 
problem involves low-dimensional dynamical system as well as the nonlinear transient dynamics 
of a jet in cross flow. Using the developed framework we compute the modes associated with 
the transient but also asymptotically unstable directions of the phase space and we assess their 
time-dependent stability properties. 
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2 Optimally time-dependent (OTD) basis for transient in¬ 
stabilities 

Let the dynamical system 


z = 

defined on a state space A C MA. We denote by St (zq) the position of the trajectory at time t 
that is initiated on zq. Also, let 

u = with L{zA) = zF{zA)^ (1) 


denote the linearized dynamical system around the trajectory St and let the inner product 
between two elements zi and ^2 be denoted as zi • Z 2 . The linear time-dependent dynamical 
system represented by equation 0 has the solution 

u{t) = (2) 

where is the propagator that maps the state of the system at time to to t. The propagator 
can be represented as the ordered product of infinitesimal propagators 




t 

to 



(3) 


where tj lies in to + {j — l)St < tj < to + jSt with t = to + n6t. Our aim is to evolve a basis 
Ui^i = 1, i.e. a set of time-dependent, orthonormal modes, so that Ui{t) optimally follows 
for all times. To achieve this goal we formulate the following quantity, which measures 
the distance between the action of infinitesimal propagator on an ortho normal basis 

Ui{t) and Ui{t + St). We have: 


F=+ St^O, (4) 

i=l 

where l/(t) = [ui (t) ,U 2 (t), ...,Ur (t)] is an arbitrary and time-dependent orthonormal basis, i.e. 
for every time instant it satisfies the orthonormality condition: 

Ui{t) ■ Uj (t) = dij, = (5) 

We observe that in the functional given by equation 0: 

Ui{t + St) = Ui{t) + Stiii + 0{St‘^), 

^t+6t _ ^LiSut)6t _ J ^ StL{StA) + 0{St^), 


where I is the identity matrix. Replacing the above equations into the functional given by 
equation Q results in: 


r 

T{ui,U 2, ...,Ur) = '^ 


dUj (t) 
dt 
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(6) 


In Figure [T] we illustrate the distance function F. For each direction Ui we have, due to the 
normalization property = 1, the rate of change Ui lying in the orthogonal complement of 
Ui. Under this constraint we choose as iii the vector that is closest to the image of Ui under the 
effect of the operate L. 

We emphasize that the minimization of the function 0 is considered only with respect to 
the time-derivative (rate of change) of the basis, t/(t), instead of the basis U{t) itself. This is 
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Figure 1: Illustration of the distance function T utilized to define the time-dependent modes. 


because we do not want to optimize the subspace that the operator is acting on, but rather find 
an optimal set of vectors, that best approximates the linearized dynamics in the subspace 

U. We then solve the resulted equations and compute U{t). We will refer to these modes as 
the optimally time-dependent (OTD) modes, and the space that these modes span as the OTD 
subspace. 

Note that a direct minimization of the function (|^, over finite-time intervals, with respect to 
the modes U{t)^ would result the well-known Euler-Lagrange equations. However, in this case 
the emphasis is put on finding the optimal ‘input’ to the operator L in order for the function ^ 
to be minimized, while here our aim is different, i.e. we seek to find an optimal set of vectors 
that approximate L in the best possible way, when L is acting on U. 

To summarize the above discussion, we are considering the following minimization problem: 

For a trajectory St find U (t) such that ^ minimized. 

In the reminder of this paper we will demonstrate that this minimization principle allows to 
capture any arbitrary and transient kind of growth caused by the linear operator, e.g. non-normal 
or exponential. In what follows we obtain a differential formulation for the evolution equations 
that correspond to the proposed minimization problem. 


2.1 Evolution equations for the time-dependent modes 

Before we proceed to the minimization of the function given by equation ([^, we express 
the orthonormality constraints for the basis elements in terms of their time derivatives. By 
differentiating over time, we have 


duj (t) 
dt 


■ Uj {t) + 


duj jt) 

dt 


Ui{t) = 0, i,j = 


The above condition is satisfied if the following condition is valid: 


(7) 


• Uj (t) = ipij (t ), ( 8 ) 

where (pij{t) is any arbitrary function for which cpij = —cpji. Clearly, for any choice of (Pij{t) the 
above condition guarantees that if {ui is initially ortho normal, it will remain orthonormal 

for all times. As we will see the choice of cpij does not change the evolved subspace. However, it 
allows for different formulations of the evolution equations. Using the constraint ^ we have the 
following Theorem. 
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Theorem 2.1 The minimization principle defined within the basis elements that satisfy the 
constraint is equivalent with the set of evolution equations 

= L{St,t)uk{t) - ^ {L{St,t)uk{t) -Uj (t) - (fkj{t))uj (t), k = (9) 

i=i 

where cpij is an arbitrary function for which cpij = —(fji. 

Proof: We first formulate the minimization problem that also takes into account the appropriate 
number of Lagrange multipliers, Xij (t) ,with i^j = 1, In this way we obtain: 

Ov if (^) i^) - V>i3 w] • 

2=1 ^ / jf = l ^ / 

( 10 ) 


We consider the first derivative with respect to each iii (t) to obtain the set of equations 

^ = 2 - L{Sut)uk{tfj +Yfkj {t)uj {t ). 

To obtain an extremum we need the right hand side of the last equation to vanish: 

= m, t)uk{t) -\y (^) ’ ( 11 ) 

i=i 

which should be solved together with condition We take the inner product of equation 0 
with mode ui {t) and obtain 

■ Ui {t) = L{St, t)uk{t) ■ Ui {t) - = ip,,i (2). 

Using the last equation and substituting Xki{t) in 0 will result in the evolution equations 
This completes the proof. 

In a more compact form the evolution equation for a finite-dimensional operator L G 

fh 

can be obtained, where we express the OTD subspace in a matrix U G whose column 

is Ui. The function cpij is correspondingly expressed in the matrix notation as G with 

^ = {Tij}i,j=i' Thus, the evolution equation for the OTD modes can be expressed as: 

pjTJ 

— =LU-U{U^LU-$), (12) 

where ( denotes the transpose of a matrix. 

Now we define the reduced operator Lr{t) G that is obtained by projecting the original 
operator onto the subspace U{t). Thus, 

Lr = U^LU. (13) 

Therefore the OTD equation could be equivalently expressed as: 

^=LU -U{Lr-^). (14) 

In what follows we will need to the definition of equivalence between two subspaces. 

Definition 2.1 The two OTD sub spaces U G and W G are equivalent at time t if 

there exists a transformation matrix R G such that U{t) = W{t)R, where R is an orthogonal 
rotation matrix, i.e. R^R = /. 
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In the following we show that the evolution of two OTD subspaces U{t) G and W{t) G 

under different choice of ^(t), which are initially equivalent (at t = 0), will remain equivalent 
for every time t > 0, i.e. U{t) = W(t)R(t)^ where R(t) G is an orthogonal rotation matrix 

governed by the matrix differential equation: 

— =R$u- ^wR, (15) 

R{0) = i?o, 


with ^jj G and G being the two different choices of ^ for the evolution of U and 

W, respectively, while Rq in the initial orthogonal rotation matrix, i.e. 1/(0) = V{0)Ro. We first 
prove the following Lemma. 


Lemma 2.1 The solution R{t) to the matrix differential equation given by equation lT3f , remains 
an orthogonal rotation matrix for every time t > 0 given that the initial eondition R{R) is an 
orthogonal rotation, i.e. i?(0)^i?(0) = I, and skew-symmetrie matrices. 


^ . d(R^R) 

Proof: We show that -;- 

dt 

d{R^R) _ 
dt 


= 0 for every t > 0 : 

R^R + R^R 

{R^u - ^wRfR + R^{R^u - ^wR) 
—^jjR^R T R^T^wR T R^R^u — R^T^wR 


= R^R^u - ^uR^R- 


Clearly Rf^R = / is a fixed point for the above equation. 

Next we prove that the for a given dynamical system, the OTD subspaces that are initially 
equivalent, remain equivalent for all times. 


Theorem 2.2 Suppose that U{t) G and W{t) G 

with different choices of <T{t) functions denoted by ^u{t) ^ 

We also assume that the two bases are initially equivalent, i.e. U{0) = W{0)Rq, where Rq G 
is an orthogonal rotation matrix. Then the subspaces U (t) and W (t) are equivalent for t > 0 
with a rotation matrix R{t) governed by the matrix differential equation (15). 


satisfy the evolution equation (12) 
and ^wif) G respectively. 

rxr 


Proof: We plug U{t) = W{t)R{t) into the OTD equation for U{t): 


WR + WR = LWR - WR{R^W'^LWR - 4>u). 

Multiplying both sides from the right with BR and using the identity RR^ = I, results in: 

W = LW - W{TW'^LW - muR'^ + RR^)- 
Now we substitute R from equation ( |l5| into the above equation to obtain: 

W = LW - W{W'^LW - $w), 

which is the evolution equation for the OTD basis W{t). This completes the proof. 

Theorem |2.2| implies that the difference between the two bases will evolve along directions 
already contained in the initially common subspace. To this end, both bases will continue to span 
the same subspace and the variation between the two is only an internal rotation. Therefore, 
the two family of equations will result in the same time-dependent subspace. There are multiple 
choices for the function cpij and we now examine a special one. 
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2.1.1 The Dynamically Orthogonal formulation 

is (fij = 0 for all i^j. The resulted evolution 


The simplest choice for the function ip in 
equations in this case will have the form 


f = 


(16) 


where Q := I — UU^ is the orthogonal projection operator onto the subspace U. Note that ipij = 0 
corresponds to the dynamical orthogonality (DO) condition |44] that has been employed to 
derive closed equations for the solution of stochastic PDEs. In this case uncertainty is resolved 
only along specific modes that evolve with time by projecting the original equation of the system 
over these directions. The evolution of these modes (stochastic subspace) is done according to 


equations derived using the DO condition and they have the general form of system (16). The 
equivalence of system (16) with the minimization problem § provides a clear interpretation for 
the evolution of the DO modes. 


2.2 Steady linearized dynamics 

Here we consider the special case where L is a time-independent operator. We prove that the basis 
defined through the introduced minimization principle will asymptotically span the eigenvectors 
of L associated with the most intense instabilities (i.e. eigenvalues with largest real part). In 
particular we have the following theorem: 

Theorem 2.3 Let L G be a steady and diagonalizable operator that represents the lin¬ 

earization of an autonomous dynamical system. Then 

i) Equation (10) has (^) = equilibrium states that consist of all the r-dimensional 

subspaees in the span of r distinet eigenveetors of L. 

a) From all the equilibrium states there is only one that is a stable solution for equation (10). 
This is given by the subspaee spanned by the eigenveetors of L assoeiated with the r eigenvalues 
having the largest real part. 

Proof: Let = d>A where is the matrix of eigenvectors of the operator L with the column of 
d> G being the eigenvectors: d> = {0i|02| • • • l^n}? and A G is the diagonal eigenvalue 

matrix whose entries are: A = diag(Ai, A 2 ,..., A^). 

i) First, we show that a subspace Uq that is in the span of precisely r eigenvectors of the 

operator L, is an equilibrium state for equation ([T^. Without loss of generality, we consider 
the first r eigenvectors to span such space, i.e. Uq £ = span{ 0 i| 02 | • • • |0r} associated 

with Ar = diag(Ai, A 2 ,..., A^), and therefore Uq can be expressed in eigenvector coordinates as: 
Uo = ^ri^o where no G is the projection coefficients. Therefore LUq = L^ri^o = ^r^ri^o ^ 

span{0i|02| • • • |0r}- As a result, LUo is in the null space of the orthogonal projector Q, and 
thus QLUo = 0. This result is independent of the choice of eigenvectors. To this end, we note 
that, for an operator L with distinct eigenvalues, there exists a number of (^) = of such 

equilibrium states. 

ii) Next, we show that from all the equilibrium states the subspace that is spanned 

by the eigenvectors associated with the largest eigenvalues, is the only stable equilibrium. 
First, we investigate the stability of Uq G We denote the complement of the space 
with = span{(/)^+i|(/)^+ 2 | • • • |0n} associated with the corresponding eigenvalues of A^c = 
diag(A^+i, A^+ 2 , • • •, An). We consider a perturbation U' e that belongs to the orthogonal 

complement of I/q, i.e. U' G l/(p: 

U{t) = UoBeU\t), U\t)TUo. 
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We note that the orthonormality condition for U{t), i.e. = /, is satisfied for e << 1: 

U{tfU{t) = (f/o + eU'{t)f{Uo + eU’{t)) 

= U^Uo + e{U^U'(t) + U'itfUo) + e^U'(tfU'(t) 

= I + e^U\tfU\t) 

/, for e << 1. 


In the above equation, we use the orthonormality condition of U' ± Uq that implies: U^U^t) = 
U'{t)'^Uo = 0. Moreover, since U^U' {t) = 0, we immediately obtain: 




T dU'jt) 
dt 


= 0 


t > 0, 


(17) 


which requires the evolution of the perturbation, i.e. dU'{t)ldt, to remain orthogonal to Uq for 
all times. Now, linearizing the evolution equation stated in equation (16) around the equilibrium 
state Uo yields: 

^ = QLU' - UoU'^LUo - U'UJLUq. 

In the above equation, the term U'^LUq = 0, since LUq G span{(/)i 1021 ... |0r} and U' ± Uq. 
Therefore the evolution equation for the perturbation equation becomes: 


BTT' 

— = QLU' - U'U^LUo. 


(18) 


Next, we transform the evolution equation (18) into the eigenvector coordinates. The perturbation 
U' can be expressed in eigenvector coordinates as: U' = where n' G Replacing 

Uq = and U' = into equation (18) results in: 


(Ik' 

$— = Q^Kk' - 


Having assumed that L is non-deficient, we multiply both sides of equation (19) by ^ 

cIk' 

= IIAk' — k' KQ^J^r-^rf^O 
= UAk' — k' KQ^ArKo^ 


(19) 


( 20 ) 


where we introduced H = ^ and used the orthonormality condition of UqUq = = 

I. We then multiply both sides by to obtain: 


dt ^ 

where p' = k'kq^. The perturbation p' can be decomposed as: 


P = 


Pr 


where p'^ and p'^c G with r + = n. The matrix H can be written as: 

n = Inxn - ^ 

~ Atxn ^ • 

In the above equation, we note that: 

Ir 


$-1$, = 


0,. 


KoKq^J {^r ^r<^) = {KoKQ^r^rKoKQ^ 

~ prxr ©rXTc) • 


( 21 ) 


( 22 ) 


(23) 


(24) 

(25) 







In the last equation, we used the orthonormality condition of UqUq = = I and 

introduced Q = kq Kq . Replacing equations (24) and (25) into equation (23) results in: 


n 


*rxr ^rXTc 
Or^xr 


By replacing equation (26) into equation (21), the evolution equation becomes: 


d 

dt 


0 0A^c\ /p; 

0 A^c j Ip 


PrA^ 


(26) 


(27) 


The evolution equation for the perturbation must satisfy the orthogonality constraint dU'{t)ldt T 
Uq expressed by equation 0- The orthogonality constraint requires that: 

= Uj{QW - U'UjLUo) = 0. 


Since Uq G the above orthogonality condition, in fact, imposes r constraints on the evolution 

of perturbation p'(t). In the following, we simplify these constraints. In the above equation: 

UAQLU' - U'U^LUo) = 

= — k'^0 

= /^o ^^4>(nAp' - p'A^)/^o 
= (4>^ 4>^c) (IIAp' — p'Ar)Ho 

= Hq (^I>pxr ^rxrc) (BAp p Arf^^HQ. 

Therefore, the orthogonality constraint dU'{t)/dt T Uo is equivalent to: 

{jrxr 0rxrc) (BAp — p A^^) = 0. (28) 


It follows that 


{jrxr 0rxrc) 


0 QAr 
0 Arc 


y]-d 


rxr ^rxr, 


Qr 


PrAr 

p'r Ar 


= 0 . 


or equivalently. 


QAr^Pr^ Pr^r — ^{prA'^ ^rcPvc)' 


(29) 


(30) 


Using the orthogonality constraint given by equation ( |3Q| ), and using equation ( |27[ ), the evolution 
equation for p(.(t) becomes: 


^=e{p',Ar-ArJJ- 


(31) 


The above equation shows that the evolution of p(.(t) can be expressed solely based on the 
evolution of p'r A)' Thus, the stability of the solution only depends on the stability of p(.^. The 
evolution of p(. in the index notation is given by: 


dp'ij 




i = r + l,r + 2 ,... ,n. 
Therefore, the asymptotic stability of p' requires that: 


real(A^) < real(Aj), for 


i = r + l,r + 2 


,..., 11 /, 


and j = 1, 2,..., r. (32) 


and j = 1, 2,..., r. (33) 


The inequality ( [^ shows that the subspace ^r with eigenvalues having the largest real part is 
the only stable solution to equation (16). This completes the proof. 
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2.3 Time-dependent linearized dynamics 


Consider the evolution of an arbitrary (not orthonormal) subspace V{t) G under the 

time-dependent linearized dynamics which is governed by: 

dV 

= L{m (34) 

4"(0) = K), 

and consider the corresponding OTD evolution: 

^ = L{t)U - ULrit), (35) 

U{0) = Uo, 

We choose the initial state of the OTD basis such that Uq and Vq span the same subspace. In 
the following Theorem we show that the OTD modes exactly span the subspace V{t). More 
specifically, we show that the two subspaces U{t) and V{t) are related via a time-dependent 
transformation matrix. 


Theorem 2.4 Let V{t) G be the subspaee evolved under the time-dependent linearized 

dynamics, eq. (^]); and U{t) G be the OTD basis evolved with eq. (35). Then, assuming 

that initially the two subspaees are equivalent, i.e. there is a matrix Tq such that Vb = UqTq, 
there exists a linear transformation T{t) such that V{t) = U{t)T{t), where T{t) is the solution of 
the matrix differential equation: 


^ = Lr{t)T{t) (36) 

r(o) = To. 

Proof: We substitute the transformation V{t) = U{t)T{t) into the evolution equation for V{t)\ 

tlT + Uf = LUT 


Substituting T from equation (36) after rearrangement results in: 

(it = LUT - ULrT. 

Multiplying both sides of the above equation by T~^ from the right results in: 

U = LU- ULr, 


which is the OTD equation for U{t). The initial condition is also U{0) = Uq. This completes the 
proof. 

For dynamical systems with persistent instabilities the evolution under the time-dependent 
linearized equation (34) is unstable, and V{t) grows rapidly. Even for stable dynamical systems, 
as we move away from t = 0, almost any vector approaches to the least stable direction. However 
the evolution of the OTD modes, due to their built-in orthonormality, is always stable, and as 
we will demonstrate in our results, the OTD evolution leads to a well-conditioned numerical 
algorithm that peels off the most unstable directions of the dynamics. 


2.4 Mode ranking within the subspace 

Having derived the subspace that corresponds to the most unstable directions, the next step is to 
rank these directions internally, i.e. within the subspace. As we describe below, there are two 
ways to rank the OTD basis based on the growth rate. Both of these approaches amount to an 
internal rotation within the OTD subspace. 
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(1) The instantaneous growth rate (7i{t) [22l[25] is obtained by computing the eigenvalues of 
the symmetric part of i.e +1/^)/2: 

(37) 

where H = diag(cri(t), cr 2 (t),..., (Tr{t))^ and R^{t) G is the matrix of eigenvectors. We rank 

these values from the least stable to the most stable directions, i.e.: 


> CF2{t) • • • > CFr{t). 


We note that crmax(0 = numerical abscissa of the operator L^. Therefore crniax(0 

represents the maximum instantaneous growth rate within the OTD subspace. 

(2) The instantaneous eigenvectors of the reduced operator can also be obtained through the 
equation: 

LrR^ = R^K^ (38) 

where the eigendirections are represented by R^{t) G The instantaneous eigenvalues are 

denoted by A^(t) = diag(Ai(t), A 2 (t),..., A^(t)), where A^ are ranked from the eigenvalue with 
the largest real part Ai(t), to the eigenvalue with the smallest real part Xr{t). 

Based on the above two rankings we define the rotated OTD basis: 

[/^’^(t) = U{t)R^’^{t), (39) 


where U^^^{t) G is the ranked representation of the OTD modes defining the space [/, and 

is the rotation matrix obtained from either of the two described strategies for mode ranking 
([ 3 ^). For a non-Hermitian operator L, the ranking based on the instantaneous 


(eq. (37) or eq. 

eigenvectors R^{t), results in modes which many not be mutually orthogonal. For this case, 
the orthogonal representation of the least stable directions can be obtained by performing a 
Gram-Schmidt orthonormalization. 


3 Linear dynamics 

Linear dynamics and in particular non-normal behaviour has a critical role in determining the 
short-time fate of a disturbance in both linear and nonlinear dynamical systems. To study the 
evolution of an initial condition under the effect of linear or linearized dynamics, reduction in 
eigenfunction coordinates is often utilized. However for highly non-normal operators, a 
large number of eigenfunctions are required to correctly capture the non-normal behavior, since 
the eigenfunctions are far from orthogonal and, in fact they constitute a highly ill-conditioned 
basis {i.e. \\U^^^\\\\{U^^^)~^\\ ^ 1). In what follows we demonstrate the computational efficiency 
of the OTD modes in capturing non-normal behavior and contrast the OTD basis with the 
eigenfunction basis for linear systems. 

3.1 Advection-diffusion operator 

As the first case, we consider the advection-diffusion operator which has a wide range of appli¬ 
cations in fluid mechanics, financial mathematics, and many other fields. Particularly we are 
interested to study the effect of non-normality on the reduced operator, both in the short-time 
and the long-term asymptotic behavior. The operator with zero boundary condition is given by: 

( 40 ) 

u{0R) = 0 , u{lR) = 0 . 

where u is the diffusion coefficient and c is the advect ion speed. The evolution equation for the 
OTD basis is expressed by: 

du ■ 

-^ = L{ui) - {L{ui),Uj)uj, i,j,= l,...,r. (41) 

Mi(0,t)=0, Ui{l,t)=0, 
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(a) (b) 

Figure 2: Instantaneous eigenvalues of the reduced operator with r = 4 for c = 1. The dashed 
lines show the eigenvalues of C with the largest real part, (a) = 0.2; (b) = 0.02. 


where the inner product and the induced norm are: 


(u^v) := / u{x)v{x)dx^ ||iz|| := 

Jo 


To solve the evolution equation we use Chebyshev collocation discretization implemented in 
chehfun [T]. The OTD basis is initialized with orthonormalized modes: 


Unix.O) = 


sin(n7rx) 

II sin(n7rx) || ’ 


n = l,2 ,...,r, 


for all the cases considered in this section. 

In the case of a nearly normal operator, i.e large z/, an optimal basis must be close to the 
dominant eigenfunctions for all times. On the other hand, for a set of parameters that corresponds 
to a non-normal operator, an optimal basis should differ from the eigenfunctions for short-time 
dynamics, and only for t ^ 1 should it converge to the operator eigenfunctions. In the remaining 
of this section, we demonstrate that the OTD subspace captures short- and long-time dynamics 
for both normal and non-normal operators. 

To analyze the behavior of the OTD basis we choose r = 4. The instantaneous eigenvalues of 
the reduced operator, real(Ai(t)), for u = 0.2 are shown with colored solid lines in Figure [^a). 
The real part of the least stable eigenvalues of the operator L are shown with the dashed lines. It 
is clear that the instantaneous eigenvalues quickly approach the least stable eigenvalues of L. In 
Figure l^a), the first two elements of the OTD basis, internally rotated in the eigendirection, 
and 02, are shown at t = 0.2 and t = 4.0. These are superimposed with the two most unstable 
eigenfunctions of the operator. At t = 0.2, the modes are very close to the eigenfunctions of L 
and at t = 4.0 the modes have essentially converged to the eigenfunctions. This demonstrates 
that due to the strongly normal behavior of the operator at u = 0.2, the eigenfunctions explain 
the dynamics accurately for both short time and long time, and the OTD basis quickly converges 
to the space spanned by the eigenfunctions. 

At z/ = 0.02, however, the instantaneous eigenvalues converge with a much slower rate and 
much later, i.e. t > 3, to the operator eigenvalues, as it is shown in Figure [^b). Accordingly, as 
it can be seen in Figure [^b), the OTD basis elements are different from the eigenfunctions of 
1/ at t = 0.2, and only until later the basis approaches the eigenfunctions. Another important 
observation is related with the advection direction of the OTD basis, which is left-ward. For 
instance at t = 0.2, the basis has advected an approximate distance of Ax = cAt = 0.2 to the 
left. As a result the OTD basis has a strong presence in the region 0 < x < 0.8. This is not the 
case for the eigenfunctions of L which are primarily concentrated near x = 0. 
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(a) 


(b) 


Figure 3: OTD basis for r = 4 and c = 1. The circles show the first two dominant eigenfunctions 
of L. The solid and dashed lines show the first and second OTD modes, i.e. u^(xA) and Un(xA): 
(a) z/ = 0.2; (b) z/ = 0.02. 


3.2 Orr-Sommerfeld/Squire operator 

As the second case we consider the Orr-Sommerfeld/Squire (OS/SQ) equation that governs 
the evolution of infinitesimal disturbances in parallel viscous flows. The eigenvalues of OS/SQ 
operator are highly sensitive to perturbations, and its eigenfunctions are linearly dependent, 
resulting in a highly ill-conditioned linear dynamical system. To this end, the OS/SQ equation 
is considered only for demonstration purposes, i.e. to illustrate that the OTD modes capture 
correctly the short-time evolution of the infinitesimal disturbances, as well as their asymptotic 
(long-term) behavior. 

We consider the plane Poiseuille flow in which the base-flow velocity is uni-directional given 
by u(x,^, z) = U{y)ex^ with U(y) = 1 — 1 /“^. The disturbance is taken to be: 


v'(x, y, z, t) = v(^, t) exp(i(ax -h i^z), 


with 


V = 


and V = 


<n J \^J 

where v' and rj' are the wall-normal velocity and the vorticity, respectively, and a and p denote 
the streamwise and spanwise wave numbers, respectively. 

The Orr-Sommerfeld equation in velocity-vorticity formulation is given by: 


dw 


with boundary conditions: 


where L is a linear operator: 


V = Vv = 77 = 0 


at 


^ = ±1, 


T — (^os 0 

\Lc LsqJ^ 


(42) 

(43) 

(44) 


with: 


Los = - k^)-^ 

Lc = -ipVU, 

Lsq = - k^) - iaVU, 


— (p2 _ ^2)2 _ ^2) 

Re 
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and k = is the modulus of the wave vector and V = d/dy. For the complete derivation 

of the OS/SQ equations, we refer to [51]. 

For the choice of the inner product we use the energy measure, which provides a physically 
motivated norm that arises naturally from the OS/SQ equation [50]. The inner product is given 
by: 

1 

(vi,V2 )^:=py v^Mv2dy, (45) 

where: 

?)' 

and ( denotes complex conjugate. In the following we consider a discrete representation of the 
operator L G Assuming a solution of the form: 


V = 0exp At, 


the initial value problem (42) transforms to an eigenvalue problem of the form 


= 4>A, 


where A = diag(Ai, A 2 ,..., A^) and 4> = {0i|02| • • • |0n} are, respectively, the matrices of the 
eigenvalues and the eigenvectors of L. 

Orszag [39] showed that for Re < Rcc — 5772.22 all eigenvalues of the operator L have 
negative real parts and therefore any perturbation is asymptotically stable. However, even for 
Re < Rec^ the energy of a perturbation may experience significant transient growth. This is a 
direct consequence of the non-normality of L. In this section we look at the evolution of the 
OTD modes for the OS/SQ operator. In particular we consider a three-dimensional perturbation 
with a = 1 and ^ = 1 at Re = 5000 that corresponds to an asymptotically stable operator. 

Since the dynamical system is linear, the linear tangent operator and L are identical. Thus, 
the evolution equation for the OTD modes U = {ui, U 2 ,..., u^} becomes: 

^ = Lui-{Lui,Uj)^Uj, i,j = 1,2,... ,r, (47) 

with the boundary conditions: 

Vi = Vvi = rji = 0 at ^ = ±1 i = 1,2,..., r. (48) 

For space we use Chebyshev collocation discretization with 256 points, while for time advancement 
we use the first order implicit Euler method. 


3.2.1 Initial condition 

We initialize the OTD subspace such that its span encompasses the optimal initial condition: an 
initial condition that reaches the maximum possible amplification at a given time t = tmax- 
The optimal initial condition can be formulated as m- 

G(^max) = max^^— j —~ || exp(Atmax) |||; (49) 

voT^o llvoll^; 

The value of || exp(Atmax)|||; is equal to the principal singular value si of the propagator 
= exp(Ltmax)- It follows that: 

^tmaxvo = VS, (50) 

where V(tmax) = {vi(tmax), Vn(tmax)} are the right singular eigenvectors and Vq = 
{vio, V 2 o,..., v^q} are the left singular eigenvectors, and S = diagjsi, S 2 ,..., consists of 
singular values of the operator B. The initial state of the subspace of size r is chosen to be: 

Uio = Vio, i = l,...,r. 


14 





The admissible initial conditions for the OTD modes must satisfy (i) the orthonormality constraint, 
and (ii) the boundary conditions at the walls given by equation ( [4^ . It is straightforward to 
show that the above choice is compatible with these criteria. We also note that short of these 
criteria, the choice of initial conditions for the OTD subspace is arbitrary. Certainly, the choice 
of optimal initial condition is of high practical importance with significant attention paid to 
in the literature. We refer the readers to m and references therein. Moreover, due to the 
non-normality of the operator, the optimal initial condition requires a large number of eigenmodes 
for accurate representation, resulting in a relatively high-dimensional system in eigenmode 
coordinates compared with the OTD modes. 


3.2.2 Transformation matrix 


We obtain a time-dependent reduction of the OS/SQ operator by projecting L onto the OTD 
subspace using the energy inner product: 


iryW = {ui,L{uj))^ i,j = 


(51) 


where Lr{t) G is the reduced OS/SQ operator. The reduced operator is then used to 

evolve the transformation matrix T(t) G according to equation (36). Using the same initial 


condition Vq for both OTD modes and OS/SQ, results in Tq being the identity matrix, i.e Tq = /. 
In the following section we compare the evolution of Vq under the full OS/SQ operator with the 
evolution of Vq using the transformation relation V(t) = U(t)T(t). 


3.2.3 Asymptotically stable subspace 

We consider the evolution of the OTD subspace with r = 2 and r = 3 for three-dimensional 
perturbations at Re = 5000 and streamwise and spanwise wave numbers of o = 1 and ^ = 1. 
At this Reynolds number all perturbations are asymptotically stable, while some perturbations 
experience significant non-normal growth in the short-time regime. In Figure the norm of the 
solution operator, ||e^^|||;, is shown. As it can be seen, the energy of some initial conditions is 
amplified by a factor of over one hundred. The maximum energy growth can be achieved at 
^max = 25.06 for the optimal initial condition. The optimal initial condition is obtained from 
equation ( [^ . We initialize the two OTD modes with the first two elements of the right singular 
eigenvectors Vq. 

Now, we first compare the evolution of the initial subspace with r = 2 using the reduced 
operator with that of the full OS/SQ operator for the choice of initial condition explained in 
Section [3.2.1 In Figure]^ we compare the energy of vi(t) and V 2 (t) obtained from (i) evolution of 
the OTD and using the transformation matrix T(t) to obtain ^i{t) = Uj(t)Tj^(t), i, j = 1,..., r, 
and (ii) by solving the full OS/SQ operator, i.e. v^(t) = ~ both cases, excellent 

agreement in both short-time and large-time evolution is observed. This demonstrates that the 
OTD modes correctly follow the evolution of a class of initial conditions according to Theorem 
2.4 Given that at Re = 5000 the OS/SQ operator is highly non-normal, a large number of 
eigenmodes are required to correctly follow the evolution of an initial condition. However the 
OTD modes, do not require additional modes beyond the dimension of the initial subspace. 
In Figure the instantaneous eigenvalues for r = 2 along with numerical abscissa for r = 2 
and r = 3 are shown. The black dashed lines show the real parts of the eigenvalues of the 
OS/SQ operator. The eigenvalues shown are the four most unstable ones of the OS/SQ operator. 
The significant non-normal energy growth manifests itself with positive real eigenvalues and 
instantaneous growth rates in finite time, despite all eigenvalues of OS/SQ having negative real 
parts. The instantaneous eigenvalues for the case with r = 2 converge to the first two least stable 
eigenvalues of the OS/SQ operator in accordance to Theorem |2.3[ For r = 2, the largest real 
instantaneous eigenvalue and the numerical abscissa amax are nearly identical. This shows that 
ui(t) is nearly aligned with the direction of maximum growth for all times. Now we explore 
some aspects of increasing the dimension of OTD from r = 2 to r = 3. For the sake of brevity, 
let us denote the quantities for the case r = 3 with the superscript ( )'. The initial condition 
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Figure 4: Transient energy growth, G{t) = ||v|||. for plane Poiseuille flow at Re = 5000 with 
a = 1 and (3 = 1. The solid lines (blue and red) show the energy growth of two different initial 
perturbations computed with the reduced operator with r = 2. The circles show the exact energy 
growth computed with the OS/SQ operator. 



Figure 5: Plane Poiseuille flow at Re = 5000, a = 1 and ^ = 1\ instantaneous real eigenvalues 
with r = 2 (blue lines); numerical abscissa with r = 2 (gray lines); numerical abscissa with r = 3 
(red); the real part of the first four least stable eigenvalues of the OS/SQ operator (dashed lines). 
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Figure 6: Snapshots of the OTD modes with r = 2 ai Re = 5000 and a = 1 and /3 = 1 in the 
streamwise plane x = n. The contour shows the vertical velocity. 


of the cases with r = 2 and r = 3 are Vq = {viq,V 2 o} and Vq = {vi^, V 2 o, vs^} respectively. 
Clearly Vq C Vq and consequently V(t) C V'(t) for all times. From the transformation between 
U(t) and V(t), we can deduce that the embedding of the initial condition is preserved for the 
OTD subspaces as well, i.e. U(t) C U'(t) for all times. Therefore it is to be expected that the 
maximum growth rate in the case of r = 3 must be always larger (or equal) than the corresponding 
value in the case of r = 2. In other words crmax(^) < <^max(^) times. This observation is 

confirmed in comparing the numerical abscissa for r = 2 and r = 3 as shown in Figure The 
abrupt changes in the values of numerical abscissa are the result of eigenvalue crossing in the 
symmetric part of the reduced operator where the direction of maximum growth switches 
from one eigendirection to the other. 

Figure shows the OTD modes for the case of r = 2 at three time instants of their evolution 
at the streamwise plane x = it: (i) the initial state t = 0, (ii) maximum energy t = tmax? and (iii) 
at large time t = 300. The initial state is marked by flow patterns that oppose the base shear. 
As times evolves from t = 0 to t < tmax, the OTD modes tilt into the mean shear flow, resulting 
in significant growth rates for the subspace. At t = 300 the modes eventually approach the two 
most unstable eigenmodes of the OS/SQ operator. This demonstrates that the time-dependent 
modes capture characteristically different regimes in the evolution of the subspace. 

4 Nonlinear dynamics 

Here we consider two nonlinear systems for which we compute the OTD modes. The first system 
is a low-dimensional dynamical system, while the second one is a more complex application 
involving an unstable flow with strongly transient characteristics. 

4.1 Low-dimensional dynamical system 

We design a low-dimensional dynamical system in order to demonstrate transient growth over 
different directions and how these can be captured by the developed approach. In particular we 


17 








































Figure 7: (a) A trajectory of the considered dynamical system colored according to the state 
variable Z 3 . The non-normal vector field for 2:3 = 0 is also shown; (b) A single cycle of the 
trajectory is shown together with a single OTD mode; (c) The time series for Zi{t)^i = 1, 2, 3; (d) 
The three eigenvalues of the linearized operator are plotted with blue dashed curves while the 
growth rate of the single OTD mode is shown with red color. 


consider the following system: 


dzi 

dt 

dZ2 

dt 

dzs 

dt 


-aiZi + 62^2 + bzs 


e ^zi - a 2 Z 2 


bzsi 




- 1 ), 


(52) 

(53) 

(54) 


where we take ai = a 2 = 2, e = 0.05, and b = 20. For these parameters the dynamical system has 
an almost periodic behavior, where each cycle exhibits four distinct regimes: ( 1 ) a non-normal 
growth in the zi — Z 2 plane, (2) exponential decay in the zi — Z 2 plane to the origin, (3) an 
exponential growth in the 2:3 direction, and (4) an exponential decay in the 2:3 direction. In 
Figure we present the trajectory of the dynamical system colored according to the variable 2:3 
in the 3D phase space. The four regimes as described above can be observed. We also present 
the projection of the vector field in the 2:3 = 0 plane, where the non-normal structure can be seen 
as well. On the other hand, the singularity at + 2:1 = 0 induces a severe exponential growth 
when the state approaches this region. This configuration allows for the repeated occurrence of 
non-normal and exponential instabilities. 

We note that due to the exponential instability close to the origin the system undergoes 
chaotic transitions between positive and negative values of 2 : 3 . In Figure we present a single 
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cycle together with a single OTD mode (red vector) shown for different time instants, while the 
time series for the state variables for one cycle is shown in Figure [^. The instantaneous growth 
rate a corresponding to the computed OTD mode is shown in Figure [^. We can clearly observe 
that the OTD mode initially captures the severe exponential growth and subsequently captures 
the non-normal growth. On the other hand the eigenvalues of the full linearized operator can only 
capture the exponential growth, even in regimes where it is not relevant, while they completely 
miss the non-normal growth. 

4.2 Vertical Jet in Crossflow 

The jet in crossflow is an important problem in fluid mechanics with a wide range of applications 
from film cooling of gas turbines, fuel injection in combustion chambers, to pollutant dispersal 
from chimneys. The interplay of the jet and crossflow creates coupled vortical structures whose 
interactions are highly unsteady and three dimensional, often leading to turbulence and resulting 
in a high dimensional dynamical system [7]. The stability of the jet in crossflow has been recently 
studied in [9], where an unstable base flow is computed by forcing the Navier-Stokes equation to 
an unstable steady solution using the selective frequency damping method [2]. The Navier-Stokes 
equation is then linearized around the base flow and the global eigenmodes of the linearized 
Navier-Stokes equation are then computed. 

In this section, we compute the OTD modes for the vertical jet in crossflow In a weakly 
turbulent regime. In particular, we follow the short- and long-time evolution of the OTD subspace 
with the time-dependent base flow^ which is obtained by performing Direct Numerical Simulation 
(DNS) of the incompressible Navier-Stokes equation. 

4.2.1 Problem specification 

The problem setup is analogous to several recent studies in the literature pig. A 2D schematic 
of a vertical jet in crossflow is shown in Figure where a vertical jet is issued into the crossflow. 
The characteristic length is the displacement thickness of the crossflow boundary layer. The origin 
of the coordinate system is placed at the center of the jet exit with the jet diameter D = 3(5*. 
The computational domain spans from x = —9.375(5* to x = 55(5* in the streamwise direction, 
from ^ = 0 to ^ = 50(5* in the wall normal direction and from 2 ; = —15(5* to 2 ; = 15(5* in the 
spanwise direction. 

At the crossflow inlet the Blasius boundary layer profile with the displacement thickness of 
6 * and free-stream velocity of Uoo is specified. The jet velocity profile is given by: 

Vjet{r) = R{1 - r^)exp(-(r/0.7)"^), (55) 

where r is the normalized distance from the center of the jet: 

r = 2/D\/x^ -h 

and R = is the ratio of peak jet velocity to the crossflow velocity. The Reynolds number, 
based on the crossflow velocity I/oo and the displacement thickness, is given by Re^ — UooS'^ ju, 
while the jet Reynolds number is given by Rcj = V^Djv. We use a velocity ratio R = 3, and a 
Reynolds number Re^o = 100 or equivalently Re^ = 900. At the top boundary the free-stream 
velocity = {t7oo5 0,0} is imposed. In the spanwise direction periodic boundary conditions are 
used. At the outflow boundary a zero-normal gradient is enforced for velocity components. 

4.2.2 OTD equations for Navier-Stokes 

To compute the time-dependent base flow, denoted by := U 5 (x, t), we solve the incompressible 
Navier-Stokes equation given by: 

^ + (Ufe ■ V)Ufe = -VP6 + ;^V2U6, 

V ■ Uft = 0, 
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Figure 8: Schematic of the vertical jet in crossflow m x — y plane. A snapshot of the base flow is 
visualized by the volume rendering of the scalar held. 



Figure 9: Instantaneous real eigenvalues and the numerical abscissa of the reduced operator 
(using the OTD modes), for the jet in crossflow at Rej = 900. 


along with the boundary conditions given in subsection |4.2.1 
modes is given by: 


The evolution equation for the 


= £Afs(Ui) - {uj,CNs{^i))^o (58) 

V • Ui = 0 


where Cns is the linearized Navier-Stokes operator, given by: 

.Cjvs(ui) = -(Ub ■ V)uj - (Ui ■ V)Ub + - Wpi. 

Re 

A zero boundary condition for u^, i = 1 ,..., r is enforced at the inflow, wall, jet exit and the top 
boundaries. Periodic boundary conditions is used in the spanwise direction and at the outflow, 
while zero-normal gradient is imposed on velocity components. In the above evolution equation 
the choice of the inner product is the second L 2 norm in the complex space: 

(u,v)= / UxVx +UyVy +UzVz, (59) 

Jn 
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where u = {ux^Uy^Uz) and v = {vx^Vy^Vz) are velocity vector fields. The reduced linear operator 
is therefore obtained from: 


= {ui,CNs{Uj)), i,j = 


(60) 


4.2.3 Initial conditions 

The initial condition for the modes is obtained by an orthonormalized space spanned by {u^(x)}^_^ 
with Ui(x) = {d^i/dy^ —0), where the two-dimensional streamfnnctions are chosen as: 

'tpi{x,y) = sm{2TTf^^x) sm{2TTfy^y)I(y), (61) 

where fx^ and fy. are the wavenumbers and T{y) is a smooth indicator function, localizing the 
modes in the main body of the jet, i.e. between ys = 1.0 and ye = 6.0. More specifically, the 
indicator function is given by: 

I{y) = tanh(( 2 / -ys)/S))- tanh(( 2 / -ye)/S)), (62) 

with S = 0.5. For the calculations that follow we choose a four-dimensional OTD basis, i.e. r = 4. 


4.2.4 Visualization 


For the visualization of the base flow we solve a passive scalar field 0 that is governed by the 
advection-diffusion equation given by: 


r)f) 

^ + (U.^V)« = 


1 

ReSc 




where Sc is the Schmidt and is chosen to be 5'c = 1 . The passive scalar is set to be ^ = 1 at 
the crossfiow inlet, 6 > = 0 at the jet inlet, periodic condition at spanwise boundaries and zero 
Neumann condition on all other boundaries. As such the jet body region is roughly determined 
by; 

jet body = {x, such that 0 < 9{x,t) < 1}. 

Moreover, by volume-rendering only selected levels of 6 >, the shear layer and vortical structures 
can be revealed. In Figure the volume rendering of 0 exposes the upper and lower shear layers 
above the jet exit, and also vortical structures further downstream. For visualizing the OTD 
modes, the iso-surface of the magnitude of velocity of the ranked OTD modes U^, colored by the 
scalar field, is shown. 


4.2.5 Numerical algorithm 

We use a spectral/hp element method to perform Direct Numerical Simulation (DNS) of the 
full Navier-Stokes equation and the evolution equation for the OTD basis. The details of the 
spectral/hp element solver (Nektar) can be found in [28]. We use an unstructured hexahedral 
mesh with 99792 elements with spectral polynomial of order four. For time integration we use the 
splitting scheme with first order explicit Euler method with time increments of At = 4 x 10“^. 
The Navier-Stokes equations are first advanced for 200 time units, by which time the nonlinear 
dynamical system has reached a statistical steady state. Due to the inherent similarities of the 
evolution equations of the OTD basis to the nonlinear Navier-Stokes equation, the computational 
cost of solving a system of r OTD modes is roughly equal to r times of a single DNS run. Since 
the base is also solved along with the OTD modes, the total computational cost is (r -h 1) times 
of a single DNS run. 


4.2.6 Non-normality and transient growth 

In Figure the instantaneous real eigenvalues and the numerical abscissa of the reduced operator 
are shown. The large disparity between the numerical abscissa and the largest real eigenvalue of 
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the reduced operator exposes a large degree of non-normality in the reduced operator The 
subspace experiences significant non-normal growth initially for 0 < t < 20. This observation is 
qualitatively in accordance with the linear stability analysis of the jet in crossflow in reference 
[9]. We refer the reader to Figure 3(c) in that article, in which the initial growth rate of the 
perturbation is much larger than its asymptotic behavior. 

The snapshots of the initial evolution of the first mode, i.e. the most unstable mode, are 
shown in Figure pT| The mode is visualized by the iso-surface of the velocity magnitude equal to 
0.03. At t = 0.4 the mode clearly captures both lower and upper shear layers. As time advances, 
the presence of the upper shear layer becomes more pronounced. This is evident in the snapshots 
in the second row in Figure pT) It should be reminded that the norm of each mode remains one 
for all times and as a result the mode quickly vanishes outside the jet body, where the magnitude 
of Ui(x, t) is significantly smaller relative to the values in the shear layer regions. 

4.2.7 Long-time behavior 

As time progresses, the subspace exhibits bursts of growth and sudden excursions into stable 
directions, as it can be seen in Figure At each time instant at least one or more unstable 
directions can be observed. These unstable directions appear either as single (real eigenvalues) or 
in pairs (complex conjugate eigenvalues). The unstable directions represent persistent instabilities 
that are the hallmark of shear flows. 

In Figure the snapshots of all four modes along with the smoke volume rendering of 
the Navier-Stokes equation are shown. The modes are visualized by iso-surfaces of the velocity 
magnitude (equal to 0.02), and colored by the scalar field 0. Each column tracks one mode at 
different time instants, starting from the top row at t = 120, with the time advancement of 
At = 2 to the next row. The modes are sorted from the most unstable directions (Mode 1) to the 
most stable directions (Mode 4). The first mode captures the vortex sheet in the upper shear 
layer of the jet. This reaffirms the strong evidence that the jet upper shear layer is unstable, 
leading to the vortex roll up further downstream [9]. The second and the third modes show 
strong presence in both the upper and lower shear layers, while the fourth mode captures the 
dominant vortical structures downstream. 

The shear layer, spanned by the OTD subspace, is a critical dynamical feature since it is 
associated with the “birth place” of the instability. The strong presence of the upper and lower 
shear layers in the large times, exposes the important role of non-normality even in the asymptotic 
dynamics of this flow. Moreover the upper shear layer remains almost steady and as such it has 
negligible turbulent kinetic energy. Therefore in energy-based reduction techniques such as Proper 
Orthogonal Decomposition (POD), the shear layer appears strongly in the time-averaged fields, 
and only weakly with modes with which unstable directions are associated. A more comprehensive 
analysis of the origin of the modes and their connection with coherent structures is currently in 
progress. 


5 Conclusions 

We have introduced a minimization formulation for the extraction of a finite-dimensional, time- 
dependent, orthonormal basis, which captures directions of the phase space associated with 
transient instabilities. The central idea is to built a set of optimally-time-dependent (OTD) 
modes with rate of change that optimally spans the vector field of the full dynamical system, 
in the neighborhood of its current state. We demonstrated how the formulated minimization 
principle can be utilized to produce evolution equations for these time-dependent modes. These 
equations require a trajectory of the system, as well as, the linearized operator and their solution 
gives a time-dependent, orthonormal basis which spans the current directions (i.e. for the current 
state of the system) associated with maximum growth. For the special case of equilibrium states 
we have shown that these modes rapidly converge to the most unstable directions of the system. 

We have demonstrated the capability of the approach on capturing instabilities caused by 
linear dynamics such as non-normal effects as well as nonlinear exchanges of energy between modes. 
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Figure 10: Initial evolution of the first mode and the trajectory of the Navier-Stokes equations, 
starting from t = 0.4 with the time advancement of At = 0.4 time units. The mode is visualized 
by the iso-surface of the velocity magnitude equal to 0.03. The time-dependent base flow (DNS) 
is visualized by smoke volume rendering of a scalar field. 


In particular we have illustrated the computation of the OTD modes in order to capture energy 
growths/exchanges occurring in: i) linear systems including the advection-diffusion operator in a 
strongly non-normal regime as well as the Orr-Sommerfeld/Squire operator, and ii) nonlinear 
systems including a low-dimensional system with both non-normal and exponential growth regimes, 
and the vertical jet in crossflow in an unstable regime. For the linear systems we demonstrated 
that the time-dependent subspace captures the strongly transient non-normal energy growth 
(in the short-time regime), while for longer times the modes capture the expected asymptotic 
behavior of the dynamics. For the low-dimensional nonlinear system we demonstrated how the 
subspace captures the most unstable directions of the dynamics, associated with exponential or 
non-normal growth, while for the fluid flow example we also explored the connection between the 
shear flow, non-normal growth and persistent instabilities. 

The proposed approach paves the way for i) the formulation of efficient, reduced-order filtering 
and prediction schemes for a variety of infinite dimensional problems involving strongly transient 
features, such as rare events, and ii) the formulation of low-energy control algorithms that will be 
able to suppress the instability in a very early stage by applying reduced-order control methods 
the moment that the instability has begun to emerge. The proposed framework should also 
be important for the fundamental understanding of the dynamical processes behind transient 
features, through the computation of finite-time Lyapunov exponents (a task that is not feasible 
in an infinite dimensional setting) and the analysis of the associated energy transfers. 
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